/*
    Copyright (C) 2009 Giacomo Spigler

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License along
    with this program; if not, write to the Free Software Foundation, Inc.,
    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*/


#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>

#include <cv.h>

#include <cutil.h>

#include "standard.h"

#include "gabor.h"




texture<unsigned char, 2, cudaReadModeNormalizedFloat> pyr_texture;  //return normalized values, as opposed to 'cudaReadModeElementType'  /  'cudaReadModeNormalizedFloat'






__global__ void S1_GPU(unsigned char *s1_0, unsigned char *s1_1, unsigned char *s1_2, unsigned char *s1_3, int *w, int *h, size_t pyr_pitch) {
  unsigned int scale, i, coord;
  int x, y, j, k;


  int w_=w[0];
  int h_=h[0];
  for(scale=0; scale<NUMSCALES; scale++) {
    int baseY=scale*h[0];
//    int w_=w[scale];
//    int h_=h[scale];
    int szI=w_*h_;

    for(i=blockIdx.x*blockDim.x+threadIdx.x; i<szI; i+=gridDim.x*blockDim.x) {
      x=i % w_;
      y=(i/w_);
      if(x>=GABOR_W && y>=GABOR_W && x<w_-GABOR_W && y<h_-GABOR_W) {

      y+=baseY;
      coord=(y-GABOR_W)*pyr_pitch+x-GABOR_W;


      float g0=0.0, g1=0.0, g2=0.0, g3=0.0, val=0.0, sum=0.0;


      for(j=0; j<GABOR_SZ; j++) {
        for(k=0; k<GABOR_SZ; k++) {
          val=tex2D(pyr_texture, x+j-GABOR_W, y+k-GABOR_W);
          sum+=val*val;
          g0+=val*gabor0[j][k];
          g1+=val*gabor45[j][k];
          g2+=val*gabor90[j][k];
          g3+=val*gabor135[j][k];
        }
      }

      sum=sqrt(sum);
      g0=abs(g0/sum);
      g1=abs(g1/sum);
      g2=abs(g2/sum);
      g3=abs(g3/sum);


//TODO: remove '*255.0' when not debugging
      s1_0[coord]=g0*255.0;
      s1_1[coord]=g1*255.0;
      s1_2[coord]=g2*255.0;
      s1_3[coord]=g3*255.0;




    }
    }


    w_/=INVERSESCALERATIO;
    h_/=INVERSESCALERATIO;
  }

}





extern "C" {


void ComputeS1_1(unsigned char *pyr_dev, unsigned char *s1_dev[4], int w_host[NUMSCALES], int h_host[NUMSCALES], int *w_dev, int *h_dev, size_t pyr_dev_pitch,  IplImage *pyr[NUMSCALES]) {

  cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<unsigned char>(); //(8, 0, 0, 0, cudaChannelFormatKindUnsigned);
/*
  pyr_texture.normalized = false;
  pyr_texture.filterMode = cudaFilterModePoint; // cudaFilterModeLinear   cudaFilterModePoint
  pyr_texture.addressMode[0] = cudaAddressModeClamp; //Wrap
  pyr_texture.addressMode[1] = cudaAddressModeClamp;*/
  cudaBindTexture2D(0, &pyr_texture, pyr_dev, &channelDesc, w_host[0], h_host[0]*NUMSCALES, pyr_dev_pitch);



  int sz=w_host[0]*h_host[0]*(NUMSCALES);
  cudaMemset(s1_dev[0], 0, sz);
  cudaMemset(s1_dev[1], 0, sz);
  cudaMemset(s1_dev[2], 0, sz);
  cudaMemset(s1_dev[3], 0, sz);


  S1_GPU<<<S1BLOCKS, S1THREADS>>>(s1_dev[0], s1_dev[1], s1_dev[2], s1_dev[3], w_dev, h_dev, pyr_dev_pitch);
  cudaThreadSynchronize();



//TO CHECK!
/*
int ori=0;
int scal=0;
cvZero(pyr[0]);
cudaMemcpy2D(pyr[0]->imageData, pyr[0]->widthStep, s1_dev[ori]+scal*pyr_dev_pitch*pyr[0]->height, pyr_dev_pitch, pyr[scal]->width, pyr[scal]->height, cudaMemcpyDeviceToHost);
*/


  cudaUnbindTexture(pyr_texture);

}






}


